#include "header.h"

static fftw_complex *a[4];
static fftw_complex *h[5];
static fftw_complex *U, *B;
static fftw_complex *bkQ[2];
static fftw_plan pf0, pb0;
static double phs1, pos1, res1, kD, kDw, kS;
static double *k, *u, overN, g, dt, nu;
static int N, reflag;
static vect rh[4], aux; 
static work_ptr tmp;
static fftw_complex *fk;
static double az;
static FILE* fhlog;
void err_msg(char* str)
{
  printf("%s\n",str);
  exit(1);
}

double signum(double x) {
  if (x > 0) return  1.;
  if (x < 0) return -1.;
  return 0.;
}

double signum2(double x) {
  if (x > 0) return  1.;
  else return -1.;
}

void init(int *argc, char **argv, params_ptr in, work_ptr wrk)
{
  char str[80], line[80], value[80], param[80];
  sprintf(str, "Usage:\n\t%s input", argv[0]);
  if (*argc != 2) err_msg(str);
  sprintf(str, "%s.cfg", argv[1]);
  FILE *fh = fopen(str,"r"); 
  fhlog = fopen("run.log","w"); 
  if (fh == NULL) err_msg("Cannot open file");
  while (fgets(line, 80, fh)!=NULL) {
    sscanf(line, "%s\t%s", param, value);
    if (strcmp(param,"#runname=") == 0) sprintf(in->runname,"%s", value);
    if (strcmp(param,"#resname=") == 0) sprintf(in->resname,"%s", value);
    if (strcmp(param,"#vrsname=") == 0) sprintf(in->resnamev,"%s", value);
    if (strcmp(param,"#npoints=") == 0) in->n = atoi(value);
    if (strcmp(param,"#kcutoff=") == 0) in->k_cut = atoi(value);
    if (strcmp(param,"#cfl_con=") == 0) in->cfl = atof(value);
    if (strcmp(param,"#gravity=") == 0) in->g = atof(value);
    if (strcmp(param,"#stokesw=") == 0) in->stokes = atoi(value);
    if (strcmp(param,"#velocit=") == 0) in->c = atof(value);
    if (strcmp(param,"#fintime=") == 0) in->tmax = atof(value);
    if (strcmp(param,"#io_time=") == 0) in->t_skip = atof(value);
    if (strcmp(param,"#skip_it=") == 0) in->skip = atoi(value);
    if (strcmp(param,"#der_cnt=") == 0) in->d_cnt = atoi(value);
    if (strcmp(param,"#readdat=") == 0) in->rflag = atoi(value);
    if (strcmp(param,"#hp_visc=") == 0) in->nu = atof(value);
    if (strcmp(param,"#readpls=") == 0) in->pflag = atoi(value);
    if (strcmp(param,"#n_poles=") == 0) in->d_poles = atoi(value);
    if (strcmp(param,"#pertres=") == 0) in->res1 = atof(value);
    if (strcmp(param,"#pertphs=") == 0) in->phs1 = atof(value);
    if (strcmp(param,"#pertpos=") == 0) in->pos1 = atof(value);
    if (strcmp(param,"#pertsft=") == 0) in->sft1 = atof(value);
    if (strcmp(param,"#vprtres=") == 0) in->res2 = atof(value);
    if (strcmp(param,"#vprtphs=") == 0) in->phs2 = atof(value);
    if (strcmp(param,"#pertpow=") == 0) in->pow = atoi(value);
  }
  fprintf(fhlog, "runname = %s\n", in->runname );
  fprintf(fhlog, "resname = %s\n", in->resname );
  fprintf(fhlog, "npoints = %d\n", in->n );
  fprintf(fhlog, "kcutoff = %d\n", in->k_cut );
  fprintf(fhlog, "cfl_con = %f\n", in->cfl );
  fprintf(fhlog, "gravity = %f\n", in->g );
  fprintf(fhlog, "stokesw = %d\n", in->stokes );
  fprintf(fhlog, "velocit = %f\n", in->c );
  fprintf(fhlog, "hp_visc = %e\n", in->nu );
  fprintf(fhlog, "t_maxim = %f\n", in->tmax );
  fprintf(fhlog, "io_time = %f\n", in->t_skip );
  fclose(fh); fclose(fhlog);
  wrk->V   = fftw_malloc(sizeof(fftw_complex)*(in->n));  
  wrk->R   = fftw_malloc(sizeof(fftw_complex)*(in->n));
  wrk->Z   = fftw_malloc(sizeof(fftw_complex)*(in->n));
  wrk->Q   = fftw_malloc(sizeof(fftw_complex)*(in->n));
  wrk->Q_t = fftw_malloc(sizeof(fftw_complex)*(in->n));
  wrk->Q_u = fftw_malloc(sizeof(fftw_complex)*(in->n));
  wrk->Z_u = fftw_malloc(sizeof(fftw_complex)*(in->n));
  wrk->Phi = fftw_malloc(sizeof(fftw_complex)*(in->n));
  wrk->P_u = fftw_malloc(sizeof(fftw_complex)*(in->n));
  for (int j = 0; j < 4; j++) {
    rh[j].r = fftw_malloc(sizeof(fftw_complex)*(in->n));
    rh[j].v = fftw_malloc(sizeof(fftw_complex)*(in->n));
  }
  U = fftw_malloc(sizeof(fftw_complex)*(in->n));
  B = fftw_malloc(sizeof(fftw_complex)*(in->n));
  bkQ[0] = fftw_malloc(sizeof(fftw_complex)*(in->n));
  bkQ[1] = fftw_malloc(sizeof(fftw_complex)*(in->n));
  aux.r = fftw_malloc(sizeof(fftw_complex)*(in->n));
  aux.v = fftw_malloc(sizeof(fftw_complex)*(in->n));
  for (int j = 0; j < 4; j++) { 
    a[j] = fftw_malloc(sizeof(fftw_complex)*(in->n));
    h[j] = fftw_malloc(sizeof(fftw_complex)*(in->n));
  }
  h[4] = fftw_malloc(sizeof(fftw_complex)*(in->n));
  u = fftw_malloc(sizeof(double)*(in->n));
  k = fftw_malloc(sizeof(double)*(in->n));
  for (int j = 0; j < in->n; j ++) {
    u[j] = 2*pi*j/(in->n);
    k[j] = j;
    if (j > (in->n)/2) k[j] = j - (in->n);
    //wrk->R[j] = 1 + ((in->res1)*cexp(1.I*pi*(in->phs1)))/( ctan(0.5*(u[j] - pi)) - 1.I*(in->pos1) ) ;
    //wrk->Z_u[j] = 1./(wrk->R[j]);
    wrk->R[j] = 1 + 0.5*(in->res1)*cexp(1.I*pi*(in->phs1))*(1./ctan(0.5*(u[j]-pi-1.I*(in->pos1))) - 1.I);
    //wrk->Z_u[j] = 1./(wrk->R[j]);
    //wrk->Z_u[j] = (ctan( 0.5*(u[j] - pi)) - 1.I*(in->res1))/(ctan( 0.5*(u[j] - pi)) - 1.I*(in->pos1)*cexp(1.I*pi*(in->phs1)));
    //wrk->R[j] = 1./(wrk->Z_u[j]);
    //wrk->Z_u[j] = 1. + 1.I*((in->res1) - (in->pos1)*cexp(1.I*pi*(in->phs1))  )/cpow(ctan(0.5*(u[j]-pi)) - 1.I*(in->pos1)*cexp(1.I*pi*(in->phs1)), in->pow);
    //wrk->Z_u[j] = 1./wrk->R[j];

    //wrk->V[j] = (in->res2)*(cexp(1.I*pi*(in->phs2)))/(tan(0.5*(u[j] - pi)) - 1.I*(in->pos1) ); // 0.05
    wrk->V[j] = 0.5*(in->res2)*cexp(1.I*pi*(in->phs2))*( 1./ctan(0.5*(u[j]-pi-1.I*(in->pos1))) - 1.0I);
  }
  pf0 = fftw_plan_dft_1d(in->n, a[0], a[1], FFTW_FORWARD, fmode);
  pb0 = fftw_plan_dft_1d(in->n, a[1], a[0], FFTW_BACKWARD, fmode);
  N = in->n; overN = 1./N; kD = 7.*(in->n)/16; kDw = 0.25*(in->n)/16; kS = 9.*(in->n)/16.;
  g = in->g; dt = (in->cfl)*2./N; 
  in->w_ptr = wrk; nu = in->nu;
  in->num_iter = floor((in->tmax)/dt) + 1;
  dt = (in->tmax)/(in->num_iter); in->dt = dt;
  phs1 = in->phs1; res1 = in->res1; pos1 = in->pos1;

  /*memcpy(a[0], wrk->R, sizeof(fftw_complex)*N);
  fftw_execute(pf0);
  printf("Zero mode of R is %e\t%e\n", overN*creal(a[1][0]), overN*cimag(a[1][0]));
  a[1][0] = 0.;
  for (int j = 0; j < N ; j++) a[1][j] = overN*a[1][j];
  fftw_execute(pb0);
  memcpy(wrk->R, a[0], sizeof(fftw_complex)*N);
  for (int j = 0; j < N ; j++) wrk->R[j] += 1.;*/

  if (in->rflag == 1) read_restart(in, wrk);

  memcpy(a[0], wrk->V, sizeof(fftw_complex)*N);
  fftw_execute(pf0);
  printf("Zero mode of V is %e\t%e\n", overN*creal(a[1][0]), overN*cimag(a[1][0]));
  a[1][0] = 0.;
  for (int j = 0; j < N ; j++) a[1][j] = overN*a[1][j];
  fftw_execute(pb0);
  memcpy(wrk->V, a[0], sizeof(fftw_complex)*N);

  reflag = 0;
  if (in->pflag == 1) read_poles(in, wrk);    
  reflag = refine(wrk, in, 0);
  if (reflag == 1) {
    printf("Check if solution is resolved ... No\n");
    err_msg("Restart with more Fourier modes\n");
  }
  printf("Check if solution is resolved ... Yes\n"); 
  update_dfunc(wrk);
  reconstruct_surface(wrk);
  fk = fftw_malloc((in->n)*sizeof(fftw_complex));
  for (int jj = 0; jj < in->n; jj++) fk[jj] = 1.; //0.5*(1+tanh((k[jj] + kD)/kDw));
}

void read_restart(params_ptr in, work_ptr wrk) {
 FILE *fhin = fopen(in->resname, "r");
 char str[256], line[256], xx[256], yy[256], dum1[256], dum2[256];
 double x, y;
 int m;

 printf("Starting from %s\n", in->resname); m = 0;
 while (fgets(line, 256, fhin)!=NULL) {
    if (m == 0) {
      fgets(line, 256, fhin);
      fgets(line, 256, fhin);
      fgets(line, 256, fhin);
    }
    sscanf(line, "%s\t%s\t%s\t%s\n", dum1, dum2, xx, yy);
    x = atof(xx); y = atof(yy);
    wrk->Z[m] = x + 1.I*y; 
    //printf("%.15e\t%.15e\n",creal(wrk->Z[m]),cimag(wrk->Z[m]));
    m++;    
 }
 der(wrk->Z, wrk->Z_u);
 fclose(fhin);
 fhin = fopen(in->resnamev, "r");
 printf("Starting from %s\n", in->resnamev); m = 0;
 if (fhin == NULL) {
   err_msg("Could not open vrsname");
 }
 while (fgets(line, 256, fhin)!=NULL) {
    if (m == 0) {
      fgets(line, 256, fhin);
      fgets(line, 256, fhin);
      fgets(line, 256, fhin);
    }
    sscanf(line, "%s\t%s\t%s\t%s\n", dum1, dum2, xx, yy);
    x = atof(xx); y = atof(yy);
    wrk->V[m] = x + 1.I*y; 
    m++;    
 }
 fclose(fhin);
 memcpy(a[0], wrk->V, sizeof(fftw_complex)*N);
 fftw_execute(pf0);
 printf("Zero mode of V is %e\t%e\n", overN*creal(a[1][0]), overN*cimag(a[1][0]));
 a[1][0] = 0.;
 for (int j = 0; j < N ; j++) a[1][j] = overN*a[1][j];
 fftw_execute(pb0);
 memcpy(wrk->V, a[0], sizeof(fftw_complex)*N);

 for (int j = 0; j < in->n; j ++) {
   wrk->Z_u[j] = 1. + wrk->Z_u[j];
   wrk->R[j] = 1./(wrk->Z_u[j]);
 }
}

void read_poles(params_ptr in, work_ptr wrk) {
  FILE *fhin = fopen(in->resname, "r");
  char str[256], line[256], xx[256], yy[256], grbg[256];
  double chi[in->d_poles], gam[in->d_poles], y0;
  char *dummy;
  
  printf("Starting from %s\n", in->resname);
  printf("Velocity set to %.15e\n", in->c);
  printf("Reading pole data ...\n");
  dummy = fgets(line, 256, fhin);
  printf("%s", line);
  dummy = fgets(line, 256, fhin);
  printf("%s", line);
  sscanf(line, "# N = %s\tL2_residual = %s\ty0 = %s\n", grbg, grbg, xx);
  y0 = strtod(xx, NULL);
  printf("y0 set to %.15e\n", y0);
  dummy = fgets(line, 256, fhin);
  
  for (int j = 0; j < in->d_poles; j++) {
    dummy = fgets(line, 256, fhin);
    sscanf(line, "%s\t%s\t%s\t%s\n", xx, yy, grbg, grbg);
    chi[j] = strtod(xx, NULL);
    gam[j] = strtod(yy, NULL);
  }
  
  for (int p = 0; p < in->n; p++) {
    wrk->Z[p] = 1.I*y0;
    wrk->Z_u[p] = 1.;
    for (int j = 0; j < in->d_poles; j++) {
      wrk->Z[p] = wrk->Z[p] + gam[j]/(tan(0.5*(u[p] - pi)) - 1.I*chi[j]);
      wrk->Z_u[p] = wrk->Z_u[p]  - 0.5*gam[j]/cpow( (tan(0.5*(u[p] - pi)) - 1.I*chi[j])*cos(0.5*(u[p] - pi)), 2);
    }
    wrk->R[p] = 1./wrk->Z_u[p];
    wrk->P_u[p] = (in->c)*(wrk->Z_u[p] - 1.);  // parameters angle pi/2, residue 1e-4, pos 0.01i
    wrk->V[p] = 1.I*(wrk->P_u[p])*(wrk->R[p]);
    wrk->Z_u[p] = wrk->Z_u[p] + 0.5*res1*cexp(1.I*pi*phs1)/cpow(sin(0.5*( u[p]-(in->sft1)*pi)) - 1.I*pos1*cos(0.5*(u[p] - (in->sft1)*pi)),2); //normal perturbation
    //wrk->Z_u[p] = wrk->Z_u[p] + 4.e-6I*(1.0e+3)*cexp(-1.I*1.0e+3*u[p]);
    wrk->R[p] = 1./wrk->Z_u[p];
  }
  memcpy(a[0], wrk->V, sizeof(fftw_complex)*N);
  fftw_execute(pf0);
  a[1][0] = 0.; 
  for (int j = 0; j < N; j++) a[1][j] = overN*a[1][j];
  fftw_execute(pb0);
}


void compute_inverse(fftw_complex *in, fftw_complex *out) {
  fftw_complex ovr_r0 = 1./in[0];
  fftw_complex incr;
  for (int j = 1; j < N/2; j++) {
    a[0][j] = out[N-j];
    a[1][j] = in[N-j];
  }
  a[0][0] = ovr_r0;
  a[1][0] = in[0];
  for (int j = 1; j < N/2; j++) {
    incr = 0;
    for (int m = 0; m < j; m++) { 
      incr += a[0][m]*a[1][j-m];
    }
    a[0][j] = -ovr_r0*incr;
  }
  for (int j = 1; j < N/2; j++) {
    out[N-j] = a[0][j];
    out[j] = 0;
  }
  out[0] = a[0][0]; 
}

void compute_zu(work_ptr in) {
  //for (int j = 0; j < N; j++) a[2][j] = in->R[j] + 1.;
  fft(in->R, a[2]);
  for (int j = 0; j < N; j++) a[2][j] = a[2][j]/N;
  compute_inverse(a[2], a[3]);
  memcpy(in->Z_u, a[3], sizeof(fftw_complex)*N);
  ifft(a[3], in->Z_u);
}

void der(fftw_complex *in, fftw_complex *out) {
  memcpy(a[0], in, sizeof(fftw_complex)*N);
  fftw_execute(pf0);
  for (int j = 0; j < N ; j++) a[1][j] = 1.I*k[j]*overN*a[1][j];
  fftw_execute(pb0);
  memcpy(out, a[0], sizeof(fftw_complex)*N);  
}

void der4(fftw_complex *in, fftw_complex *out) {
  memcpy(a[0], in, sizeof(fftw_complex)*N);
  fftw_execute(pf0);
  for (int j = 0; j < N ; j++) a[1][j] = cpow(1.I*k[j],12)*overN*a[1][j];
  fftw_execute(pb0);
  memcpy(out, a[0], sizeof(fftw_complex)*N);  
}

void hilbert(fftw_complex *in, fftw_complex *out) {
  memcpy(a[0], in, sizeof(fftw_complex)*N);
  fftw_execute(pf0);
  for (int j = 0; j < N ; j++) a[1][j] = 1.I*signum(k[j])*overN*a[1][j];
  fftw_execute(pb0);
  memcpy(out, a[0], sizeof(fftw_complex)*N);  
}

void project(fftw_complex *in, fftw_complex *out) {
  memcpy(a[0], in, sizeof(fftw_complex)*N);
  fftw_execute(pf0);
  for (int j = 0; j < N ; j++) a[1][j] = 0.5*(1.-signum(k[j]))*overN*a[1][j];
  fftw_execute(pb0);
  memcpy(out, a[0], sizeof(fftw_complex)*N);  
}
void ffilter(fftw_complex *in) {
  memcpy(a[0], in, sizeof(fftw_complex)*N);
  fftw_execute(pf0);
  for (int j = 0; j < N ; j++) a[1][j] = 0.5*(1.-signum2(k[j]))*overN*a[1][j]*cexp(-dt*overN*nu*pow(k[j],12));
  fftw_execute(pb0);
  memcpy(in, a[0], sizeof(fftw_complex)*N);  
}

void dealias(fftw_complex *in) {
  memcpy(a[0], in, sizeof(fftw_complex)*N);
  fftw_execute(pf0);
  for (int j = 0; j < N ; j++) a[1][j] = 0.5*(1.-signum2(k[j]))*overN*a[1][j]*fk[j];
  fftw_execute(pb0);
  memcpy(in, a[0], sizeof(fftw_complex)*N);  
}

void pproject(fftw_complex *in, fftw_complex *out) {
  memcpy(a[0], in, sizeof(fftw_complex)*N);
  fftw_execute(pf0);
  for (int j = 0; j < N ; j++) a[1][j] = 0.5I*k[j]*(1.-signum(k[j]))*overN*a[1][j];
  fftw_execute(pb0);
  memcpy(out, a[0], sizeof(fftw_complex)*N);  
}

void project_combined(fftw_complex *in, fftw_complex *out1, fftw_complex *out2) {
  memcpy(a[0], in, sizeof(fftw_complex)*N);
  fftw_execute(pf0);
  for (int j = 0; j < N ; j++) {
    a[1][j] = 0.5*(1.-signum(k[j]))*overN*a[1][j];
    a[2][j] = 1.I*k[j]*a[1][j];
  }
  fftw_execute(pb0);  
  fftw_execute(pb1);


  memcpy(out1, a[0], sizeof(fftw_complex)*N);
}


double find_y0(work_ptr in) {
  for (int j = 0; j < N; j++) a[0][j] = cimag(in->Z_u[j]);
  fftw_execute(pf0);
  for (int j = 0; j < N; j++) a[1][j] = -1.I*a[1][j]*overN/(k[j]);
  a[1][0] = 0; az = 0; 
  for (int j = 0; j < N; j++) az += -fabs(k[j])*pow(cabs(a[1][j]),2);
  return az;
}

void update_dfunc(work_ptr in) {
  //memcpy(bkQ[0], bkQ[1], N*sizeof(fftw_complex));
  //memcpy(bkQ[1], in->Q, N*sizeof(fftw_complex));
  compute_zu(in);
  for (int j = 0; j < N; j++) {
    //in->Z_u[j] = 1./(in->R[j]);
    in->Q[j] = 1.i*(in->V[j])/(in->R[j] - 1.);
    in->Q_t[j] = 3.*in->Q[j] - 4.*bkQ[1][j] + bkQ[0][j];
    in->Q_t[j] = in->Q_t[j]/(2*dt);
  }
  der(in->Q, in->Q_u);
  reconstruct_surface(in);
}

void reconstruct_surface(work_ptr in) {
  in->y0 = find_y0(in);
  memcpy(a[0], in->Z_u, N*sizeof(fftw_complex));
  fftw_execute(pf0);
  for (int j = 0; j < N; j++) a[1][j] = -1.I*a[1][j]*overN/k[j];
  a[1][0] = 0 + 1.I*(in->y0);
  fftw_execute(pb0);
  memcpy(in->Z, a[0], N*sizeof(fftw_complex));
  for (int j = N; j > -1; --j) {
    in->Z[j] = in->Z[j] - creal(in->Z[0]); 
  }
}

void reconstruct_potential(work_ptr in) {
  for (int j = 0; j < N; j++) a[0][j] = -1.I*(in->V[j])*(in->Z_u[j]);
  memcpy(in->P_u, a[0], N*sizeof(fftw_complex));
  fftw_execute(pf0);
  for (int j = 0; j < N; j++) a[1][j] = -1.I*a[1][j]*overN/k[j];
  a[1][0] = 0.;
  fftw_execute(pb0);
  memcpy(in->Phi, a[0], N*sizeof(fftw_complex));
}


void fft(fftw_complex *in, fftw_complex *out) {
  memcpy(a[0], in, N*sizeof(fftw_complex));
  fftw_execute(pf0);
  memcpy(out, a[1], N*sizeof(fftw_complex));
}

void ifft(fftw_complex *in, fftw_complex *out) {
  memcpy(a[1], in, N*sizeof(fftw_complex));
  fftw_execute(pb0);
  memcpy(out, a[0], N*sizeof(fftw_complex));
}

void prep_control_params(params_ptr in, work_ptr wrk) {
  for (int j = 0; j < N; j++) {
    wrk->Z_u[j] = 1./(wrk->R[j]);
  }
  reconstruct_potential(in->w_ptr);
  reconstruct_surface(in->w_ptr);
  mean_level(in);
  momentumx(in);
  momentumy(in);
  hamiltonian(in);
}

void mean_level(params_ptr in) {
  work_ptr tmp = in->w_ptr;
  double du = 2.*pi/N;
  in->current.M = 0;
  for (int j = 0; j < N; j ++) {
    in->current.M += du*cimag(tmp->Z[j])*creal(tmp->Z_u[j]); 
  }
}

void momentumx(params_ptr in) {
  work_ptr tmp = in->w_ptr;
  double du = 2.*pi/N;
  in->current.Px = 0;
  for (int j = 0; j < N; j ++) {
    in->current.Px += du*creal(tmp->Phi[j])*creal(tmp->Z_u[j]); 
  }
}

void momentumy(params_ptr in) {
  work_ptr tmp = in->w_ptr;
  double du = 2.*pi/N;
  in->current.Py = 0;
  for (int j = 0; j < N; j ++) {
    in->current.Py += du*cimag(tmp->Phi[j])*cimag(tmp->Z_u[j]); 
  }
}

void hamiltonian(params_ptr in) {
  work_ptr tmp = in->w_ptr;
  double du = 2.*pi/N;
  in->current.H = 0;
  for (int j = 0; j < N; j ++) {
    in->current.H += -0.5*du*(creal(tmp->Phi[j])*cimag(tmp->P_u[j])); 
    in->current.H += 0.5*g*du*(pow(cimag(tmp->Z[j]), 2)-0.*pow(tmp->y0, 2))*creal(tmp->Z_u[j]);
  }
}

void rhs(work_ptr in, int m) {
  memcpy(h[0], in->R, N*sizeof(fftw_complex));
  memcpy(h[1], in->V, N*sizeof(fftw_complex));
  for (int j = 0; j < N; j++) {
    h[2][j] = 2.*creal(h[1][j]*conj(h[0][j]));
    h[3][j] = h[1][j]*conj(h[1][j]);
  } 
  der(in->R, h[0]);
  project (h[2], U);
  pproject(h[2], h[1]);
  pproject(h[3], B);
  der(in->V, h[2]);

  for (int j = 0; j < N; j++) {
    rh[m].r[j] = 1.I*(U[j]*h[0][j] - h[1][j]*(in->R[j])) - 0.*nu*h[4][j];  
    rh[m].v[j] = 1.I*(U[j]*h[2][j] - B[j]*(in->R[j])) + g*(in->R[j] - 1.) - 0.*nu*h[3][j];      
  }
}

void rhs_st(work_ptr in, int m) {
  //memcpy(h[0], in->R, N*sizeof(fftw_complex));
  //memcpy(h[1], in->V, N*sizeof(fftw_complex));
  
  for (int j = 0; j < N; j++) {
    Q[j] = csqrt(R[j]);
    U[j] = 2.*creal((in->V[j])*conj(in->R[j]));
    
  }
  project(U, U);
  der(Q, h[0]);
  der(U, h[3]);
  der(in->R, h[1]);
  der(in->V, h[2]);
  for (int j = 0; j < N; j++)    
    B[j] = creal((in->V[j])*conj(in->V[j]) + 4.*(in->sigma)*cimag(h[0][j]*conj(Q[j])));
  }
  pproject(B,B);

  for (int j = 0; j < N; j++) {
    rh[m].r[j] = 1.I*(U[j]*h[1][j] - h[1][j]*(in->R[j]));
    rh[m].v[j] = 1.I*(U[j]*h[2][j] - B[j]*(in->R[j])) + g*(in->R[j] - 1.);
  }
}


void rk4(params_ptr in) {
  tmp = in->w_ptr;
  memcpy(aux.r, tmp->R, N*sizeof(fftw_complex));
  memcpy(aux.v, tmp->V, N*sizeof(fftw_complex)); 
  rhs(tmp,0);
  for (int j = 0; j < N; j++) {
    tmp->R[j] = aux.r[j] + 0.5*dt*rh[0].r[j];
    tmp->V[j] = aux.v[j] + 0.5*dt*rh[0].v[j];
  }
  rhs(tmp, 1);
  for (int j = 0; j < N; j++) {
    tmp->R[j] = aux.r[j] + 0.5*dt*rh[1].r[j];
    tmp->V[j] = aux.v[j] + 0.5*dt*rh[1].v[j];
  }
  rhs(tmp, 2);
  for (int j = 0; j < N; j++) {
    tmp->R[j] = aux.r[j] + dt*rh[2].r[j];
    tmp->V[j] = aux.v[j] + dt*rh[2].v[j];
  }
  rhs(tmp, 3);
  for (int j = 0; j < N; j++) {
    tmp->R[j] = aux.r[j] + (1./6)*dt*(rh[0].r[j] + 2.*rh[1].r[j] + 2.*rh[2].r[j] + rh[3].r[j]);
    tmp->V[j] = aux.v[j] + (1./6)*dt*(rh[0].v[j] + 2.*rh[1].v[j] + 2.*rh[2].v[j] + rh[3].v[j]);
  }
  ffilter(tmp->R);
  ffilter(tmp->V);

  memcpy(bkQ[0], bkQ[1], N*sizeof(fftw_complex));
  memcpy(bkQ[1], tmp->Q, N*sizeof(fftw_complex));
  for (int j = 0; j < N; j++) {
    tmp->Q[j] = 1.i*(tmp->V[j])/(tmp->R[j] - 1.);
  }

}

void generate_derivatives(work_ptr in, int m) {
  char filename[80];
  sprintf(filename, "der_%02d.dat", 0);
  reconstruct_surface(in);
  write_out(in->Z, filename, 0.);
  for (int j = 1; j < m + 1; j ++) {
    der(in->Z, in->Z);
    sprintf(filename, "der_%02d.dat", j);
    write_out(in->Z, filename, 0.);
  }
}

int array_max (fftw_complex *input_array, int N) {
	int i, max_i;
	double max_el;
	max_el = 0;
	max_i = 0;
	for (i = 0; i<N; i++){
		if (creal(input_array[i]*conj(input_array[i]))>max_el){
		max_el = creal(input_array[i]*conj(input_array[i]));
		max_i = i;
		}
	}
	return max_i;
}

fftw_complex* expand_complex (fftw_complex *iarray, int Nnew) {
  fftw_free(iarray);
  fftw_complex *darray = fftw_malloc (Nnew*sizeof(fftw_complex));
  memset (darray, 0, Nnew*sizeof(fftw_complex));
  
  return darray;
}

double* array_enlarge_double (double *input_array, int Nold, int Nnew) {
	double *doubled_array;

	doubled_array = (double *) fftw_malloc (Nnew*sizeof(double));
        memset ((void*) doubled_array, 0, Nnew*sizeof(double));
	fftw_free ((void*) input_array);

	return doubled_array;
}

void expand_static(int Nnew) {
  fftw_free(u); u = fftw_malloc(sizeof(double)*Nnew);
  fftw_free(k); k = fftw_malloc(sizeof(double)*Nnew);

  fftw_free(bkQ[0]);  bkQ[0] = fftw_malloc(sizeof(fftw_complex)*Nnew);
  fftw_free(bkQ[1]);  bkQ[1] = fftw_malloc(sizeof(fftw_complex)*Nnew);
  for (int j = 0; j < Nnew; j ++) {
    u[j] = 2*pi*j/Nnew;
    k[j] = j;
    if (j > Nnew/2) k[j] = j - Nnew;
  }
  fftw_free(U);  
  fftw_free(B);
  U = fftw_malloc(Nnew*sizeof(fftw_complex));
  B = fftw_malloc(Nnew*sizeof(fftw_complex));   
  memset(U, 0, Nnew*sizeof(fftw_complex));
  memset(B, 0, Nnew*sizeof(fftw_complex));
  for (int j = 0; j < 5; j++) {
    fftw_free(h[j]);
    h[j] = fftw_malloc(Nnew*sizeof(fftw_complex));
    memset(h[j], 0, Nnew*sizeof(fftw_complex));
  }
  for (int j = 0; j < 4; j++) {
    fftw_free(rh[j].r);
    fftw_free(rh[j].v);
    rh[j].r = fftw_malloc(Nnew*sizeof(fftw_complex));
    rh[j].v = fftw_malloc(Nnew*sizeof(fftw_complex));
    memset(rh[j].r, 0, Nnew*sizeof(fftw_complex));
    memset(rh[j].v, 0, Nnew*sizeof(fftw_complex));
  }
  fftw_free(aux.r);  
  fftw_free(aux.v);  
  fftw_free(fk);
  fk = fftw_malloc(Nnew*sizeof(fftw_complex));
  aux.r = fftw_malloc(Nnew*sizeof(fftw_complex));
  aux.v = fftw_malloc(Nnew*sizeof(fftw_complex));
  memset(aux.r, 0, Nnew*sizeof(fftw_complex));
  memset(aux.v, 0, Nnew*sizeof(fftw_complex));
  memset(fk, 0, Nnew*sizeof(fftw_complex));
  kD = 7.*Nnew/16; kDw = 0.25*Nnew/16; kS = 3.*Nnew/4; kS = 9.*Nnew/16.;
  for (int jj = 0; jj < Nnew; jj++) fk[jj] = 1.; //0.5*(1+tanh((k[jj] + kD)/kDw));
}

void expand_wrk(work_ptr wrk, int Nnew) {
   fftw_free(wrk->R);   fftw_free(wrk->V);   fftw_free(wrk->Q_t);
   fftw_free(wrk->Z);   fftw_free(wrk->Q);   fftw_free(wrk->Z_u);
   fftw_free(wrk->Phi); fftw_free(wrk->P_u); fftw_free(wrk->Q_u);

   wrk->R = fftw_malloc (Nnew*sizeof(fftw_complex));
   wrk->V = fftw_malloc (Nnew*sizeof(fftw_complex));
   wrk->Z = fftw_malloc (Nnew*sizeof(fftw_complex));
   wrk->Q = fftw_malloc (Nnew*sizeof(fftw_complex));
   wrk->Q_t = fftw_malloc (Nnew*sizeof(fftw_complex));
   wrk->Q_u = fftw_malloc (Nnew*sizeof(fftw_complex));
   wrk->Z_u = fftw_malloc (Nnew*sizeof(fftw_complex));
   wrk->P_u = fftw_malloc (Nnew*sizeof(fftw_complex));
   wrk->Phi = fftw_malloc (Nnew*sizeof(fftw_complex));
   memset(wrk->R, 0, Nnew*sizeof(fftw_complex));
   memset(wrk->V, 0, Nnew*sizeof(fftw_complex));
   memset(wrk->Z, 0, Nnew*sizeof(fftw_complex));
   memset(wrk->Q, 0, Nnew*sizeof(fftw_complex));
   memset(wrk->Q_t, 0, Nnew*sizeof(fftw_complex));
   memset(wrk->Q_u, 0, Nnew*sizeof(fftw_complex));
   memset(wrk->Z_u, 0, Nnew*sizeof(fftw_complex));
   memset(wrk->P_u, 0, Nnew*sizeof(fftw_complex));
   memset(wrk->Phi, 0, Nnew*sizeof(fftw_complex));
}

int refine(work *wrk, params *in, int j){
   int max_i, Nold, window;
   char s1[80];
   double value, rq, rq2, rq3, vq;
   double value2;
   Nold = in->n;

   fft(wrk->R, a[2]);       
   fft(wrk->V, a[3]);
   rq2 = 0.; rq3 = 0.; window = 5;
   for (int j = window; j > -window; j--) {
     rq2 += creal(a[2][(int)floor(kS)+j]*conj(a[2][(int)floor(kS)+j])); 
     rq3 += creal(a[3][(int)floor(kS)+j]*conj(a[3][(int)floor(kS)+j])); 
   }
   rq2 = rq2/(4*window*window);
   rq3 = rq3/(4*window*window);
   rq = sqrt(rq2/creal(a[2][0]*conj(a[2][0])));
   vq = sqrt(rq3/creal(a[2][N-1]*conj(a[2][N-1])));
   //rq2 = sqrt(creal(a[2][(int)floor(kS)]*conj(a[2][(int)floor(kS)]))/(1e-32+creal(a[2][N-1]*conj(a[2][N-1]))));
   if (j % 100 == 0) {
     fhlog = fopen("run.log","a");
     fprintf(fhlog, "Log of max/min ratio for R: %f and for V: %f at K-sample = %d\n", -log10(rq), -log10(vq), (in->n)-(int)floor(kS) );
     //printf("Log of max/min ratio for R: %f and for V: %f at K-sample = %d\n", -log10(rq), -log10(vq), (in->n)-(int)floor(kS) );
     fclose(fhlog);
   }
   if (rq !=0.0) {
      value = -log10(rq);
      value2 = -log10(vq);
      if ((value < 12.5)||(value2 < 12.5)) { 
            
      in->n = 2*Nold; 
      fhlog = fopen("run.log","a");
      //printf("Log of max/min ratio for R: %f and for V: %f\nExpanding Fourier Spectrum to %d modes\nRefining\n", value, value2, in->n);
      fprintf(fhlog, "Log of max/min ratio for R: %f and for V: %f\nExpanding Fourier Spectrum to %d modes\nRefining\n", value, value2, in->n);
      fclose(fhlog);

      expand_wrk(wrk, in->n);
      expand_static(in->n);
      a[0] = expand_complex(a[0], in->n);
      a[1] = expand_complex(a[1], in->n);
      fftw_plan pi0 = fftw_plan_dft_1d(in->n, a[0], a[0], FFTW_BACKWARD, fmode);
      fftw_plan pi1 = fftw_plan_dft_1d(in->n, a[1], a[1], FFTW_BACKWARD, fmode);
      for (int l = Nold-1; l > Nold/2; l--) {
         a[1][Nold+l] = a[2][l]/Nold;
         a[0][Nold+l] = a[3][l]/Nold;
      }
      a[1][0] = a[2][0]/Nold;
      a[0][0] = a[3][0]/Nold;

      fftw_execute(pi0); fftw_destroy_plan(pi0);
      fftw_execute(pi1); fftw_destroy_plan(pi1);
      for (int l = 0; l < in->n; l++) {
         wrk->V[l] = a[0][l];
         wrk->R[l] = a[1][l];
      }
      free_prep(in); init_prep(in);
      a[2] = expand_complex(a[2], in->n);
      a[3] = expand_complex(a[3], in->n);
      fftw_destroy_plan(pf0);
      fftw_destroy_plan(pb0);
      pf0 = fftw_plan_dft_1d(in->n, a[0], a[1], FFTW_FORWARD, fmode);
      pb0 = fftw_plan_dft_1d(in->n, a[1], a[0], FFTW_BACKWARD, fmode);
      N = in->n; overN = 1./N;

      in->dt = 0.5*(in->dt); dt = in->dt;
      in->num_iter = j + floor(((in->tmax)-in->current.T)/in->dt) + 0;
      return 1;
      } else return 0;
   } else return 0;
}


















